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1.  Introduction 

The  real  purpose  of  epidemic  theory  is  not  to  develop  interesting  and  elegant 
mathematics,  though  this  may  be  a  delightful  incidental  byproduct,  but  is  to 
facilitate  the  practical  prevention  or  control  of  actual  outbreaks  of  serious 
contagious  disease.  This  purpose  is  still  a  long  way  from  being  achieved  to  any 
appreciable  extent.  The  developed  countries  are  today  free  from  disasters  of  the 
magnitude  of  the  Black  Death  in  the  14th  century  when  perhaps  as  much  as 
25  per  cent  of  the  population  in  Europe  perished.  Nevertheless,  widespread 
epidemics  on  a  massive  scale  are  still  common  in  Africa  and  the  Far  East.  As 
the  volume  and  speed  of  modern  travel  continue  to  increase  there  is  an  ever 
growing  risk  of  the  transmission  of  virulent  infections  to  regions  where  natural 
immunity  may  be  low  though  public  health  control  is,  for  ordinary  purposes, 
more  or  less  adequate.  Even  within  a  developed  country  there  are  possible 
dangers  from  such  factors  as  the  appearance  of  new  strains  of  infectious  organ¬ 
isms  resistant  to  standard  drugs  and  antibiotics,  or  increases  in  the  contact 
rate  between  individuals  due  to  greater  population  densities  or  changes  in 
social  behavior.  The  current  increase  in  venereal  infections  in  many  countries 
could  be  a  case  in  point.  It  follows  therefore  that  it  is  eminently  worth  consider¬ 
ing  in  what  directions  research  should  proceed  in  order  to  have  an  improved 
chance  of  attaining  its  object. 

As  with  applications  to  many  other  fields  in  biology  and  medicine,  the  attempt 
to  develop  mathematical  theories  of  epidemics  exhibits  the  usual  conflict  between 
insight  and  realism.  Epidemic  theory  falls  into  two  distinct,  though  complemen¬ 
tary,  parts.  On  the  one  hand,  there  is  the  study  of  small  groups  like  individual 
families.  From  these  it  is  possible,  though  not  often  done,  to  collect  detailed 
data  to  which  relatively  realistic  models  can  be  fitted,  yielding  information 
about  such  biological  or  clinical  entities  as  contact  rate,  length  of  latent  and 
infectious  periods  and  so  forth.  However,  little  can  be  deduced  from  this  about 
the  spread  of  infection  through  a  community.  The  latter  requires  the  special 
theory  of  large  groups.  This  provides  some  insight  into  the  behavior  of  popula¬ 
tion  outbreaks,  with  regard  to  such  features  as  threshold  phenomena  or  the 
general  graphical  appearance  of  “epidemic  eurves.,,  Unfortunately,  even  the 
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simplest  mathematical  models,  especially  if  they  are  stochastic  formulations, 
present  formidable  problems  of  analysis.  Any  insights  obtained  therefore  are 
liable  to  be  most  true  of  highly  oversimplified  models.  As  soon  as  an  appreciable 
degree  of  realism  is  introduced,  especially  when  geographical  spread  is  involved, 
the  treatment  becomes  largely  intractable.  The  simplified  mathematical  theory 
is  worth  pursuing  because  of  the  hints  and  leads  it  may  produce,  but  its  limita¬ 
tions  should  not  be  forgotten. 

It  is  within  this  context  that  more  serious  consideration  should  be  given  to 
the  potentialities  of  simulating  epidemic  processes,  particularly  as  the  power  and 
availability  of  automatic  electronic  computers  are  constantly  increasing. 

Most  real  communities  entail  the  spatial  distribution  of  individuals  as  an 
essential  ingredient,  with  susceptibles  in  the  “neighborhood”  of  an  infective 
being  more  likely  to  contract  his  disease  than  those  that  are  distant.  “Neighbor¬ 
hoods”  may,  however,  be  very  complex  and  hard  to  define,  while  various  forms 
of  complicated  social  and  geographical  stratification  may  also  exist.  With  the 
aid  of  modern  computer  methods  of  data  processing  and  numerical  analysis  it 
seems  just  within  the  bounds  of  possibility  that  a  highly  diversified  community 
could  be  represented  by  a  moderately  realistic  model  and  subjected  to  inves¬ 
tigations  of  a  simulation  of  Monte  Carlo  type. 

We  shall  look  briefly  at  existing  theory,  first  for  deterministic  spatial  models, 
and  then  for  the  corresponding  stochastic  versions.  Next,  we  shall  examine 
some  new  simulation  studies  of  stochastic  epidemics  in  two  dimensions,  using 
models  which  so  far  are  highly  oversimplified,  but  are  easily  capable  of  consider¬ 
able  modification  and  extension.  Lastly,  we  shall  discuss  the  general  implications 
of  this  type  of  work  to  the  practical  problems  facing  public  health  authorities. 

2.  Deterministic  theory 

It  is  always  worth  examining  a  deterministic  model  first,  even  when  wre  are 
sure  that  a  stochastic  formulation  is  essential  to  adequate  realism,  provided 
that  the  inherent  limitations  are  not  lost  sight  of.  Stochastic  models  are  in¬ 
variably  harder  to  handle,  and  this  is  specially  true  of  epidemic  processes. 
Approximate  treatments  may  be  facilitated  if  we  know  in  advance  what  kind 
of  features  to  look  for,  and  important  indications  may  be  provided  by  an  analysis 
of  the  corresponding  deterministic  process. 

There  are  various  ways  in  which  spatial  or  topographical  elements  can  be 
incorporated  in  an  epidemic  model.  Consider,  for  example,  the  following  formu¬ 
lation  due  to  D.  G.  Kendall  (see  discussion  in  Bartlett  [3]).  We  assume  the 
existence  of  an  infinite  twro  dimensional  continuum  of  population  with  a  density 
of  a  individuals  per  unit  area.  Take  any  point  P  surrounded  by  a  small  areal 
element  dS.  Let  the  numbers  of  individuals  in  this  area  who  are  susceptible, 
infectious  or  removed  be  axdS,  crydS  and  azdS,  respectively,  where  the  quan¬ 
tities  x,  y  and  2  are  proportions  which  sum  to  unity  though  they  may  be  func¬ 
tions  of  time  and  position.  Let  the  contact  rate  be  (5  and  the  removal  rate  be  7. 
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Now  the  basic  differential  equations  defining  the  process  can  be  written  in 
the  form 

dx 
Tt 


=  —  (3  a  xy, 


(2.1) 


dy  _ 

—  =  (3axy  -  7/7, 

dz 

at  =  yv ’ 


where  y  is  a  spatially  weighted  average  of  y  given  by 

(2.2)  y(P,  t )  =  fJ\(PQ)y(Q,  t )  dS, 

in  which  dS  is  an  areal  element  at  Q  and  \(PQ)  is  an  appropriate  nonnegative 
weighting  coefficient.  Equations  (2.1)  are  an  obvious  spatial  extension  of  the 
usual  deterministic  general  epidemic  equations  (for  example,  Bailey  [1]  equa¬ 
tions  (4.5)),  in  which  the  rate  of  new  infections  is  made  to  depend  on  y  rather 
than  y. 

A  suitable  initial  condition  in  the  present  case  is  to  assume  a  focus  of  infection 
uniformly  spread  over  some  small  circle  centered  at  the  origin.  An  appreciable 
amount  of  mathematical  discussion  is  possible,  though  the  development  is  not 
without  difficulty.  So  far  as  any  practical  reference  is  concerned  the  main 
result  available  is  the  following.  Subject  to  certain  conditions,  it  can  be  shown 
that  a  pandemic  affecting  every  part  of  the  plane  will  ensue  if  and  only  if  the 
population  density  of  susceptibles  a  exceeds  the  familiar  threshold  p  =  y/(3. 
Moreover,  if  there  is  a  pandemic,  its  severity  f  is  the  unique  positive  root  of 

(2.3)  f  =  1  -  exp  (-<tf/p), 

meaning  that  the  fraction  of  individuals  eventually  contracting  the  disease  will 
be  at  least  f  in  any  area  of  the  plane,  no  matter  how  far  from  the  initial  focus. 
We  thus  have  a  pandemic  threshold  theorem  corresponding  to  the  well  known 
nonspatial  version  of  Kermack  and  McKendrick  [7]. 

An  alternative  model,  that  is  in  some  ways  more  specific,  has  been  used  by 
Bartlett  [2]  to  discuss  the  behavior  of  recurrent  epidemics.  Bartlett’s  treatment 
is  fairly  general  and  involves  terms  representing  the  migration  of  both  suscep¬ 
tibles  and  infectives,  as  well  as  the  appearance  of  new  infectives.  All  three  of 
these  features  can  be  eliminated  for  the  purpose  of  the  present  discussion,  and 
the  following  simple  argument  results.  Writing  x  and  y  for  the  actual  spatial 
densities  of  susceptibles  and  infectives,  we  assume  that  the  action  of  infection 
is  local  and  isotropic.  This  allows  us  to  write  the  first  two  differential  equations 
in  the  form 


=  ~(3x(y  +  aV-y), 


dx 
~dt 

=  fW.v  +  -  v> 


(2.4) 
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where  x  =  #(£,  y,  t),  y  =  ?/(£,  y,  t )  and  £,  y  are  the  spatial  coordinates  them¬ 
selves,  while  V2  =  d2/d£2  -j-  d2/dy2. 

In  the  initial  stages  of  an  epidemic  x  will  be  approximately  constant.  The 
second  equation  in  (2.4)  can  then  be  written  as 

(2.5)  m=Ay  +  BV2y’ 

where  A  =  fix  —  y,  B  =  afix.  Equation  (2.5)  is  of  course  a  standard  diffusion 
equation  with  solution 

(2'6)  »  =  2§t“p(^~lW!)’ 

the  constant  C  being  determined  by  the  initial  conditions. 

If  we  consider  the  total  volume  of  infection  outside  a  circle  of  radius  R ,  this 
may  be  calculated  as 

(2.7)  yR  =  j^^ydidr, 

=  2tC  exp  ( At  -  j|t)- 

Thus, 

(2.8)  R  =  2 (AB)'iH  [l  - 

As  t  — >  =o  so  R  — >  2 (AB)ll2t.  From  this  it  follows  that  the  circle  of  radius  R  cor¬ 
responding  to  any  arbitrary  value  of  ijr  grows,  for  sufficiently  large  t,  at  the 
constant  rate  R/t  =  2 (AB)112.  This  can  be  interpreted  as  the  velocity  of  'propaga¬ 
tion  of  the  disease  from  the  initial  focus. 

So  far  as  practical  consequences  are  concerned,  Kendall’s  model  raises  the 
important  issue  of  pandemic  thresholds,  and  is  in  fact  a  first  step  forward 
towards  a  more  exact  understanding  of  this  type  of  phenomenon.  Bartlett’s 
model,  on  the  other  hand,  allows  one  to  look  more  closely  at  the  actual  rate  of 
spread  of  infection.  The  approximation  derived  above  is,  however,  only  of 
limited  application  since  we  have  assumed  that  the  epidemic  is  still  in  its  early 
stages  with  roughly  constant  x. 

More  recently,  Kendall  [6]  investigated  the  deterministic  theory  of  epidemic 
spread  for  linear  communities  in  one  spatial  dimension,  and  discussed  the  forms 
of  epidemic  waves  traveling  out  from  a  focus.  The  existence  of  such  waves 
required  the  density  of  susceptibles  to  be  above  a  certain  threshold  value. 
Similar  results  might  be  expected  for  the  two  dimensional  case,  but  an  exact 
analysis  is  not  yet  available. 

No  doubt  the  more  thorough  examination  of  deterministic  models  will  in  due 
course  reveal  further  important  properties.  The  main  point  of  such  work  is, 
however,  to  see  what  light  may  be  shed  on  the  behavior  of  more  realistic  sto¬ 
chastic  versions. 
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3.  Stochastic  theory 

No  general  theoretical  treatment  is  yet  available  for  any  stochastic  analogue 
of  the  two  dimensional  deterministic  models  described  in  the  previous  section, 
though  Morgan  and  Welsh  [8]  have  recently  discussed  a  rather  special  kind  of 
stochastic  infection  process  on  a  lattice.  Bartlett  [2]  has  obtained  an  appropriate 
partial  differential  equation  for  the  probability  generating  functional  of  a  suit¬ 
able  “point”  process,  but  this  has  so  far  proved  intractable  to  a  general  analysis. 
Some  progress  was,  however,  possible  with  the  initial  stages  of  an  epidemic 
when  the  numbers  of  susceptibles  could  be  regarded  as  approximately  constant. 
Bartlett  showed  that  in  this  special  case  the  behavior  of  the  average  number  of 
infectives  at  any  time  entailed  the  propagation  of  a  wave  of  infection  similar 
to  that  found  in  the  deterministic  situation.  No  information  has  been  obtained 
on  the  general  form  of  the  epidemic  spread  as  the  stock  of  susceptibles  becomes 
appreciably  depleted. 

More  recently  Neyman  and  Scott  [9]  have  obtained  some  very  interesting 
results  for  a  two  dimensional  model  that  incorporates  a  number  of  realistic 
features  previously  neglected.  In  particular,  the  number  of  susceptibles  infected 
by  an  infective  is  made  to  depend  on  the  latter’s  location,  and  also  an  individual 
infected  at  any  point  is  allowed  to  move  away  and  become  infectious  elsewhere. 
On  the  other  hand,  it  is  a  salient  feature  of  epidemic  processes  that  each  new 
infection  diminishes  the  number  of  susceptibles  at  risk.  Neyman  and  Scott’s 
model  does  not  include  this  aspect,  and  though  approximately  valid  at  the  start 
of  an  epidemic  would  become  progressively  less  so  as  the  outbreak  built  up. 
How  serious  the  limitation  is  for  the  conclusions  reached  is  hard  to  say.  For 
example,  one  result,  under  fairly  broad  conditions,  is  that  the  probability  of 
an  epidemic  building  up  in  a  small  community  is  the  same  as  the  probability 
of  an  explosive  outbreak  over  the  whole  of  a  much  larger  area.  (Compare 
Kendall’s  pandemic  threshold  theorem  described  in  section  2.)  This  result  cor¬ 
responds  with  the  usual  public  health  view  that  neglect  of  sanitary  conditions 
in  any  part  of  a  country  may  expose  the  whole  country  to  danger.  One  might 
expect  the  limitation  referred  to  above  to  have  less  effect  on  this  conclusion 
since  in  a  nonspatial  general  epidemic  we  can  approximate  the  early  stages  by 
a  birth  and  death  process  for  the  population  of  infectives,  ignoring  the  depletion 
of  susceptibles.  And  the  probability  of  only  a  small  outbreak  is  roughly  the 
probability  of  this  process  suffering  extinction.  Another  result  that  needs  more 
careful  consideration  in  relation  to  certain  assumptions  made  concerns  the  effect 
of  an  immunization  campaign  on  the  total  size  of  epidemic.  Subject  to  certain 
conditions,  it  appears  that  the  immunization  of  a  random  proportion  0  of  the 
population  would  reduce  the  expected  size  of  epidemic  to  a  value  less  than 
(1  —  0)/0.  This  result  certainly  seems  optimistic  (as  remarked  by  Neyman  and 
Scott),  since  it  implies  that  if  only  ten  per  cent  of  the  population  were  immunized 
the  average  outbreak  would  be  less  than  nine,  a  circumstance  that  looks  im¬ 
probable  if  a  disease  were  highly  contagious  with  a  small  relative  removal  rate. 

As  a  complementary  alternative  to  theoretical  studies,  simulation  investiga- 


242 


FIFTH  BERKELEY  SYMPOSIUM:  BAILEY 


tions  can  be  undertaken.  Little  has  been  done  so  far  for  models  incorporating 
spatial  elements.  One  Monte  Carlo  study  by  Bartlett  [3]  used  a  6  X  6  grid  of 
cells.  Standard  nonspatial  continuous  time  or  discrete  time  processes  were 
assumed  within  each  cell,  while  a  stochastic  movement  of  infectives  between 
cells  with  a  common  boundary  was  adopted.  Thus,  the  spatial  element  was 
introduced  to  a  certain  extent,  though  the  model  was  primarily  concerned  with 
the  behavior  of  recurrent  epidemics  in  situations  where  there  wras  a  constant 
accession  of  new  susceptibles. 

It  was  therefore  thought  worthwhile  carrying  out  a  simulation  study  of  a 
population  of  individuals  all  of  whom  were  spread  out  spatially.  This  is  de¬ 
scribed  in  the  following  section. 

4.  Simulated  epidemics  in  two  dimensions 

In  this  section  we  describe  one  of  the  simplest  possible  models  for  a  two 
dimensional  epidemic  and  present  a  number  of  results  obtained  by  straight¬ 
forward  Monte  Carlo  simulations  carried  out  on  an  electronic  computer.  As 
emphasized  again  later,  although  this  model  is  very  much  oversimplified,  con¬ 
siderable  modifications  and  extensions  could  be  made  with  only  relatively  minor 
changes  in  the  computer  program.  Precisely  what  further  work  is  worth  doing 
is  a  matter  for  discussion. 

Let  us  envisage  a  square  lattice  of  points  each  occupied  by  a  susceptible 
individual.  It  is  convenient  to  regard  this  community  as  having  finite  size.  For 
the  present  investigation  a  square  boundary  was  chosen,  centered  on  the  origin, 
given  by  the  lines  x  =  ±  k,  y  =  ±  k.  The  total  population  size  n  is  thus  n  = 
(2 k  +  l)2.  Most  calculations  were  performed  with  k  =  5,  though  some  were 
done  for  n  =  10.  Large  populations  would  need  to  be  explored  on  a  bigger 
computer  than  the  one  actually  used,  an  Elliott  803. 

In  the  present  study  it  was  supposed  that  an  epidemic  was  initiated  by  the 
susceptible  at  the  origin  becoming  infected.  This  involves  a  certain  amount  of 
symmetry,  which  is  convenient  in  a  first  investigation,  but  not  of  course  essen¬ 
tial.  We  could  thus  examine  the  way  in  which  the  epidemic  spread  out  from  the 
focus  without  edge  effects  due  to  the  boundaries  occurring  in  the  initial  stages, 
and  ensuring  that  their  influence  was  symmetrical  when  they  did  occur. 

A  discrete  time  model  of  chain  binomial  type  (see  Bailey  [1]  for  general 
discussion  and  references)  was  adopted  involving  two  forms,  one  representing 
a  simple  epidemic  with  no  recovery  and  one  a  general  epidemic  with  infectives 
undergoing  removal.  We  first  describe  the  simple  epidemic. 

4.1.  Simple  epidemic  until  no  removal.  Let  us  suppose,  as  in  ordinary  chain 
binomial  theory,  that  after  infection  there  is  a  fixed  latent  period  following 
which  the  infective  becomes  highly  infectious  for  a  very  short  period  of  time. 
During  this  infectious  phase,  ideally  contracted  to  a  point,  there  is  probability 
p  of  any  susceptible  at  risk  contracting  the  disease  from  the  infective  in  question. 
But  a  given  susceptible  may  be  exposed  to  several  infectives.  If  the  latter  are  r 
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ill  number,  the  probability  of  the  susceptible  becoming  infected  is  1  —  (1  —  p)r. 
In  the  chain  binomial  theory  of  Greenwood  or  Reed-Frost  type  an  infective  is 
assumed  to  recover,  or  at  least  be  removed,  immediately  after  the  point  of 
infectiousness.  But  since  in  the  simple  epidemic  we  have  no  removal,  a  conven¬ 
ient  assumption  is  that  the  infective  again  becomes  highly  infectious  after  a 
further  interval  equal  to  the  original  latent  period,  and  so  on  indefinitely.  In 
short,  each  infective  has  a  series  of  infectious  points,  following  the  initial  infec¬ 
tion,  all  separated  by  intervals  equal  to  the  latent  period.  The  epidemic  will 
therefore  spread  in  a  series  of  discrete  stages  or  generations. 

Now  we  also  have  to  decide  what  individuals  are  potentially  at  risk  from  any 
given  infective.  As  we  have  deliberately  introduced  a  spatial  element  by  spread¬ 
ing  out  the  population  over  a  lattice  some  restriction  is  evidently  needed,  other¬ 
wise  we  should  simply  be  returning  to  a  nonspatial  model.  The  most  obvious 
assumption  would  be  to  regard  only  the  four  nearest  neighbors  as  being  at  risk 
(unless  already  infected).  This  type  of  assumption  is  often  used  in  various  kinds 
of  random  walk  models,  and  might  well  turn  out  to  be  mathematically  trac¬ 
table.  However,  there  are  certain  disadvantages.  One  is  that  such  a  restriction 
might  be  thought  unduly  strong,  though  this  is  perhaps  a  small  point  in  a 
preliminary  study.  More  serious  is  the  limitations  involved  in  using  a  small 
computer. 

It  is  easy  to  see  that  if  the  disease  is  so  infectious  that  p  =  1,  then  it  will 
spread  in  a  deterministic  manner  covering  at  the  {/th  generation  a  complete 
square  whose  corners  are  at  the  four  points  (±g,  ±g).  As  this  square  has  its 
diagonals  along  the  main  x  and  y  axes,  the  edge  effects  which  start  occurring 
when  g  —  k  +  1  are  more  elaborate  than  they  would  be  if  the  square  had  its 
sides  parallel  to  the  main  axes  so  as  to  be  similarly  situated  to  the  population 
boundaries.  Of  course  we  could  change  the  latter  by  turning  the  population 
square  through  45°,  but  this  makes  the  computer  program  much  more  compli¬ 
cated  if  we  are  fully  to  utilize  a  given  square  array,  the  elements  of  which  repre¬ 
sent  the  positions  of  the  individuals  in  the  population. 

Accordingly,  it  seemed  simpler  to  assume  that  the  eight  nearest  neighbors  to 
any  infective  were  at  risk  (unless  already  infected),  as  shown  in  figure  1.  With 
this  arrangement,  the  deterministic  spread  occurring  when  p  —  1  means  that 
at  the  gth  generation  a  complete  square  is  infected  whose  sides  are  given  by  the 
lines  x  =  ±<7,  y  =  ±g.  Thus  at  generation  k  the  population  square  is  com¬ 
pletely  covered,  whereas  at  generation  k  —  1  none  of  the  susceptibles  on  the 
boundary  is  infected.  The  edge  effects  are  now  very  simple  in  form.  The  epi¬ 
demic  simply  spreads  steadily  from  the  focus  until  the  boundary  is  reached, 
when  it  ceases. 

Of  course  we  are  not  specially  interested  in  the  case  p  =  1,  but  for  smaller 
values  of  p  we  may  expect  a  probability  spread  of  disease  in  which  some  degree 
of  symmetry  is  retained  vis-a-vis  the  boundaries. 

The  computed  simulations  were  programmed  in  ALGOL  and  carried  out  on 
the  small  Elliott  803  computer  in  the  Unit  of  Biometry,  Oxford.  The  general 
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scheme  of  the  computations  may  be  briefly  outlined  as  follows.  An  array  {an} , 
with  i,  j  =  —  k,  —k  4-  1,  •  •  ■  ,0,  •  •  •  ,  k  —  1 ,  k,  is  used  to  represent  the  popula¬ 
tion  of  individuals.  If  the  individual  at  the  point  (i,j)  is  susceptible  o,-y  =  0, 
if  infectious  =  1.  And,  initially,  g  =  0  with  a0 0  =  1,  all  other  a,y  being 
zero.  At  the  gth  generation  all  parts  of  the  array  that  could  have  been  reached 
by  the  epidemic,  that  is,  anywhere  in,  or  on  the  boundaries  of,  the  square  formed 
by  the  lines  x  =  db g,  y  =  dag,  are  systematically  inspected.  If  any  a,j  =  1,  no 
change  is  possible.  But  if  a,j  =  0,  we  have  a  susceptible  who  may  be  at  risk. 
All  the  eight  nearest  positions  must  be  examined  and  the  sum  of  the  values  of 


Figure  1 

Individuals  at  risk  from  a  given  infection. 


the  relevant  an  formed.  Let  this  be  r.  If  r  =  0  there  are  no  adjacent  infect ives. 
If  r  >  0  then  a,j  must  be  changed  from  0  to  1  with  probability  1  —  (1  —  p)r. 
An  appropriate  pseudorandom  number  is  calculated,  and  the  relevant  transition 
performed.  The  random  number  algorithm  xn+i  =  ox„(mod  235)  based  on 
Behrenz  [4]  was  used.  This  is  very  quick  to  compute  though  not  entirely  satis¬ 
factory  from  a  statistical  point  of  view  because  of  appreciable  serial  correlations 
between  successive  terms.  It  was  however  considered  to  be  sufficiently  accurate 
for  the  purpose  in  hand.  We  proceed  in  this  way  for  some  arbitrary,  but  suffi¬ 
ciently  large  number  of  generations,  in  order  to  complete  a  single  simulation. 

In  the  investigations  whose  results  are  described  below  25  separate  simula¬ 
tions  were  performed  for  each  epidemic  set  up  with  a  different  value  of  p.  While 
a  larger  number  of  repetitions  is  desirable,  say  100,  the  results  seemed  generally 
useful  with  many  coefficients  of  variation  being  of  the  order  of  10  per  cent  or 
less.  It  would  be  very  easy  to  run  much  larger  numbers  of  simulations  on  a  big 
computer  using  the  same  ALGOL  programs.  Also,  a  more  satisfactory  random 
number  generator  could  be  inserted  in  the  program. 

Now  it  would  be  tedious  to  interpret  and  wasteful  of  computer  storage 
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capacity,  to  try  to  keep  records  of  the  stochastic  behavior  of  the  process  at  each 
point  of  the  lattice.  We  may  conjecture  that  some  of  the  symmetry  present 
when  p  =  1  still  remains  for  smaller  values  of  p.  In  which  case  it  is  convenient 
to  amalgamate  results  for  the  boundaries  of  a  square  of  given  size.  There  is  bound 
to  be  some  lack  of  symmetry  at  the  corners  since  the  risk  of  infection  is  always 
rather  less  there,  but  we  have  chosen  to  ignore  this.  When  any  stage  of  a  partic¬ 
ular  simulation  has  been  completed,  we  record  the  total  sum  of  the  changes  in 
the  an  values  along  the  boundaries  of  each  square,  and  also  keep  a  running- 
total  of  these  sums  and  a  running  total  of  the  squares  of  these  sums  for  all 
simulations  to  date.  At  the  end  of  the  set  of  25  simulations  we  can  calculate 
means  and  standard  errors  for  the  quantity  in  question.  This  quantity  is  in  fact 
an  epidemic  curve,  on  the  usual  definition  (see  Bailey  [1]),  for  the  whole  set  of 
individuals  along  the  boundaries  of  a  given  square. 

Results  are  shown  for  the  case  lc  =  5,  that  is,  an  11  X  11  square,  in  figures 
2,  3,  4,  and  5,  with  p  =  0.2,  0.4,  0.6  and  0.8,  respectively.  The  figures  shown 
against  the  graphs  themselves  refer  to  the  five  available  squares.  The  upper¬ 
most  graphs  (a)  in  each  figure  record  the  average  epidemic  curve  for  each 
square  as  a  whole,  while  the  middle  graphs  (b)  show  the  average  per  individual 
in  each  square.  This  latter  quantity  is  therefore  the  epidemic  curve  relating  to 
the  average  position  in  a  given  square.  Finally,  the  lower  graphs  (c)  give  the 
epidemic  curve  for  the  whole  population  of  121  individuals,  and  the  graphs  (d) 
show  the  distribution  of  the  completion  time,  that  is,  the  number  of  generations 
it  takes  for  the  disease  first  to  infect  everyone.  The  material  for  a  given  value 
of  p  was  produced  during  a  single  set  of  computations,  each  involving  25  rep¬ 
etitions,  averaging  about  three  hours  computing  time  in  all. 

Although  it  might  have  been  statistically  more  satisfactory  to  exhibit  these 
results  in  the  form  of  histograms,  it  would  then  have  been  impossible  to  super¬ 
impose  the  several  curves  without  confusion.  The  frequency  polygon  method 
was  therefore  preferred. 

The  (a)  curves  in  each  figure  show  how  the  epidemic  spreads  more  rapidly 
as  p  increases,  and  how  the  magnitude  of  the  epidemic’s  effect  is  greater  at 
greater  distances  from  the  focus,  though  of  course  more  time  is  required  for  a 
build  up  the  further  out  we  go.  When  p  is  small  the  curves  are  relatively  flatter 
and  the  build  up  takes  longer  than  when  p  is  large. 

Of  more  immediate  interest  are  perhaps  the  (b)  graphs  which,  as  already 
mentioned,  show  in  effect  the  epidemic  curves  relating  to  an  average  individual 
on  a  given  square.  The  curves  can  also  be  regarded  as  frequency  distributions 
since  the  areas  under  them  must  all  be  unity.  These  distributions  are  therefore 
more  directly  comparable  with  one  another  and  constitute  one  way  of  represent¬ 
ing  the  advancing  epidemic  wave.  It  will  be  seen  that  the  wave  appears,  in 
general  terms,  to  progress  at  a  more  or  less  steady  rate,  though  its  effect  is  more 
spread  as  we  move  further  from  the  focus,  at  least  for  small  p.  When  p  is  large, 
for  example,  0.8,  the  waves  are  almost  identical  in  form  irrespective  of  the 
distance  from  the  origin.  And,  of  course,  when  p  =  1  the  spread  is  entirely 
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deterministic  with  exactly  one  new  case  appearing  per  point  at  precisely  the  gth 
generation  for  all  points  on  the  square  x  =  ±g,  y  =  ±g. 

The  more  or  less  linear  spread  of  disease  can  also  be  seen  in  the  lower  (c) 


graphs  which  show  epidemic  curves  for  the  whole  population.  At  least,  the 
average  growth  is  linear  until  about  the  sixth  generation  when  edge  effects 
begin  to  appear.  As  p  increases  the  epidemic  curve  drops  more  sharply.  When 
p  =  0.2  the  fall  is  slower  than  the  rise,  but  with  p  =  0.8  the  cutoff  becomes 
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quite  steep,  due  to  the  epidemic’s  reaching  the  outer  boundary  with  a  higher 
degree  of  probability.  Again,  with  p  =  1  the  cutoff  would  be  completely  vertical. 
The  reason  for  the  linear  rise  is  presumably  that  in  the  rather  restricted  model 
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Figure  3 

Simple  epidemic  on  an  11  X  11  square  with  p  =  0.4. 


adopted  each  infective  is  in  contact  with  only  a  small  number  (eight)  of  other 
individuals,  and  there  is  not  the  same  opportunity  for  the  probability  of  infec¬ 
tion  to  rise  to  the  values  obtained  when  there  is  full  homogeneous  mixing. 
It  is  not  clear  whether  anything  useful  can  be  deduced  from  the  completion 
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Simple  epidemic  on  an  11  X  11  square  with  p  =  0.0. 


time  graphs  at  (d).  Since  the  whole  of  each  curve  is  based  on  only  25  observa¬ 
tions,  sampling  fluctuation  is  large.  But  it  can  be  seen  that  completion  time 
becomes  shorter  and  less  variable  as  the  disease  becomes  more  infectious. 

A  number  of  simulations  were  also  run  for  a  21  X  21  square  with  k  =  10. 
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Figure  5 

Simple  epidemic  on  an  11  X  11  square  with  p  =  0.8. 

The  general  conclusions  are  rather  similar  to  those  above,  though  we  can  of 
course  follow  the  progress  of  the  epidemic  for  a  greater  number  of  generations 
before  edge  effects  predominate.  As  a  single  example,  let  us  consider  the  epi¬ 
demic  curve  for  p  =  0.6  shown  numerically  in  table  I.  It  is  evident  that  for  the 
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TABLE  I 


Epidemic  Curve  for  a  Simple  Epidemic  in  a 
Whole  Population  Spread  Oct  Over  a 
21  X  21  Lattice  with  p  =  0.6 


Generation 

Average  No.  of  New  Cases 

1 

4.4 

2 

11. S 

3 

18.6 

4 

26.7 

5 

33.1 

6 

41.2 

7 

48.5 

8 

56.0 

9 

63.5 

10 

72.2 

11 

40.1 

12 

18.1 

14 

0.6 

1.) 

0.2 

first  ten  generations  there  is  an  almost  constant  increase  of  about  7.5  new  cases 
per  generation,  after  which  the  epidemic  curve  falls  away  rather  sharply  as 
before. 

4.2.  General  epidemic  with  removal.  The  simple  epidemic  model  of  the  pre¬ 
vious  section  can  readily  be  extended  to  cover  the  more  general  case  when  re¬ 
moval  of  some  kind  is  envisaged.  We  merely  have  to  adopt  the  familiar  chain 
binomial  assumption  that  after  the  first  occurrence  of  infectiousness,  the  infec¬ 
tive  in  question  ceases  to  be  infectious  and  enters  a  third  state,  represented  in 
the  simulations  by  a,j  =  2,  say.  Alterations  required  to  the  previous  computer 
program  are  minimal.  The  set  of  results  corresponding  to  the  simple  epidemics 
described  above  are  shown  in  figures  6  to  9  for  the  same  range  of  values  of  p, 
though  the  distributions  of  completion  time  are  now  omitted  and  the  distribu¬ 
tions  of  total  epidemic  size  are  considered  in  the  separate  table  II.  The  total 
computing  time  for  each  value  of  p  was  about  six  hours. 

For  p  =  0.4,  0.0,  and  0.8  it  has  been  possible  to  use  the  same  scales  as  for 
the  simple  epidemic,  but  for  p  =  0.2  large  changes  were  necessary.  This  corre¬ 
sponds  to  some  kind  of  a  threshold  between  p  =  0.2  and  p  =  0.4. 

Comparing  figure  0(a)  with  figure  2(a)  for  instance,  we  see  how  in  the  general 
case  the  curves  present  quite  a  different  appearance.  The  average  number  of 
new  cases  per  square  now  falls  off  with  distance  instead  of  conversely,  and  all  num¬ 
bers  are  absolutely  much  less  than  in  the  simple  epidemic.  It  is  clear  that  we  are 
now  dealing  with  a  small  local  outbreak  only,  whose  presence  is  felt  less  and  less 
as  we  go  farther  from  the  focus.  The  point  is  made  even  more  obvious  by  com¬ 
paring  the  “point”  epidemic  curves  of  figure  0(b)  with  figure  2(b).  In  the  general 
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epidemic  we  have  waves  traveling  from  the  center  that  are  rapidly  damped 
out.  Similarly,  figure  6(c)  shows  the  existence  of  quite  a  small  outbreak  com- 


General  epidemic  on  an  11  X  11  square  with  p  =  0.4. 

pared  with  figure  2(c)  when  we  take  into  account  a  tenfold  difference  in  the 
vertical  scales. 

When  we  come  to  look  at  figure  7  and  figure  3,  we  find  that  for  p  =  0.4  the 
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Figure  8 


General  epidemic  on  an  11  X  11  square  with  p  =  0.6. 
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Figure  9 


General  epidemic  on  an  1 1  X  1 1  square  with  p  =  0.8. 
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differences  are  much  less  marked.  The  general  epidemic  is  now  definitely  spread¬ 
ing  in  an  explosive  manner,  though  somewhat  more  slowly  and  diffusely  than  in 
the  simple  epidemic  with  no  recovery.  For  larger  values  of  p  the  differences 
between  the  general  and  simple  types  of  epidemic  become  progressively  less 
marked  so  far  as  epidemic  curves  and  epidemic  waves  of  advance  are  concerned. 

The  range  between  p  =  0.2  and  p  =  0.4  obviously  deserves  closer  attention, 
and  a  series  of  simulations  were  carried  out  at  intervals  of  0.04  with  special 
reference  to  the  distribution  of  the  total  size  of  epidemic.  These  distributions, 
fairly  coarsely  grouped,  are  shown  in  table  II.  Although  the  samples  are  all 


TABLE  II 


Distributions  of  Total  Epidemic  Size  for  Different  Values  of  p 


Total  No.  of  New  Cases 

.20 

Values  of  p 
.24  .28 

.32 

.36 

0-  20 

13 

10 

8 

2 

0 

21-  40 

6 

4 

1 

1 

0 

41-  60 

5 

6 

3 

2 

0 

61-  80 

1 

2 

8 

4 

1 

81-100 

0 

3 

5 

8 

3 

101-120 

0 

0 

0 

8 

21 

Total 

25 

25 

25 

25 

25 

rather  small,  the  steady  shift  in  distribution  with  changes  in  p  is  immediately 
apparent.  When  p  =  0.20,  just  over  half  the  epidemic  sizes  were  20  or  less  in 
a  total  population  of  121.  At  p  =  0.24  there  are  strong  indications  of  a  spread 
over  the  whole  range,  while  by  p  =  0.32  the  concentration  is  building  up  at  the 
top  end.  At  p  =  0.30,  over  80  per  cent  of  epidemics  have  a  total  size  of  more 
than  100.  It  would  be  interesting  to  examine  the  precise  form  of  this  change 
over  from  one  extreme  distribution  to  another,  using  a  much  larger  number  of 
simulations  for  greater  accuracy. 

5.  The  role  of  simulation 

Considerable  emphasis  was  placed  in  the  introduction  on  the  importance  oi 
theoretical  epidemic  studies  having  some  practical  relevance.  We  saw  in  section 
2  how  some  preliminary  insights  could  be  obtained  into  the  behavior  of  epi¬ 
demics  incorporating  spatial  elements  by  means  of  deterministic  models.  But 
when  we  examined  more  realistic  stochastic  models  in  section  3  it  seemed  that 
these  formulations  were  either  intractable,  or  could  be  managed  only  when 
confined  to  the  opening  stages  of  an  epidemic,  or  could  be  developed  only  if 
fairly  severe  simplifications  were  assumed.  Section  4  on  the  other  hand  dealt 
with  simulation  of  some  very  simple  spatial  models,  in  which  the  limitations 
were  of  a  different  kind.  It  is  highly  pertinent,  therefore,  to  ask  in  what  way 
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such  simulation  studies  can  be  regarded  as  having  practical  applications,  and 
as  providing  useful  results  that  cannot  be  obtained  by  purely  mathematical 
manipulations  of  the  relevant  models. 

As  mentioned  at  the  beginning  of  section  4,  although  the  models  subsequently 
discussed  were  considerably  oversimplified  they  were  readily  capable  of  exten¬ 
sive  development  and  modification  by  means  of  comparatively  small  changes  in 
the  computer  programs  used.  Thus  we  assumed  that  individuals  were  perma¬ 
nently  attached  to  the  points  of  a  square  lattice,  though  able  to  infect  any  of 
their  eight  nearest  neighbors.  This  might  be  interpreted  as  permitting  some 
degree  of  movement,  at  least  sufficient  to  result  in  a  degree  of  contact  with 
nearest  neighbors.  There  is  no  reason  why  we  should  not  introduce  the  pos¬ 
sibility  of  more  distant  individuals  being  infected,  the  chance  of  infection 
diminishing  perhaps  with  distance  according  to  some  assumed  law.  This  means 
that,  when  at  any  stage  we  are  inspecting  the  array  {ai}}  and  find  that  a  partic¬ 
ular  dij  =  0,  that  is,  the  point  (i,  j)  is  occupied  by  a  susceptible,  we  should  have 
to  search  a  larger  area  around  the  point  (t,  j)  for  possible  infectives  than  en¬ 
visaged  in  the  simple  study  described.  More  computer  time  would  be  required 
but  no  new  principle  is  involved.  Indeed,  if  required,  all  kinds  of  heterogeneous 
spatial  structures  could  be  built  into  the  model.  Moreover,  additional  stages 
might  be  introduced  into  the  fundamental  stochastic  process  so  as  to  relax  the 
assumptions  of  constant  latent  period  and  point  infectiousness. 

How  far  can  one  go  in  developing  extensions  before  even  the  largest  existing 
computer  is  inadequate  is  difficult  to  say.  Perhaps  a  more  important  question 
at  the  present  time  is  what  would  be  practically  useful  if  the  computations 
could  be  performed  in  reasonable  periods  of  time.  There  are  two  obvious  pos¬ 
sibilities. 

First,  far  more  highly  structured  models  could  be  developed  that  contained 
representations  of  a  very  large  number  of  realistic  features.  The  probable  con¬ 
sequences  of  a  variety  of  public  health  measures  such  as  immunization  cam¬ 
paigns,  the  use  of  quarantining  and  restriction  of  public  movements,  and  so 
forth,  could  then  be  ascertained  by  reference  to  models  that,  though  hypothet¬ 
ical,  were  realistic  and  typical  in  the  sense  of  incorporating  a  large  number  of 
actual  clinical,  biological,  social  or  geographical  features.  The  general  insights 
obtained  from  such  work  might  be  very  much  more  powerful  than  the  elementary 
notions  of  thresholds  at  present  available.  They  might  also  usefully  be  combined 
with  the  application  of  operational  gaming  methods  to  public  health  training- 
programs.  These  have  been  employed,  for  instance,  by  the  California  Depart¬ 
ment  of  Public  Health  together  with  the  Systems  Development  Corporat  ion  [5], 
and  entail  the  construction  of  a  simulated  city  called  “Epiville.” 

Secondly,  it  might  be  possible  to  construct  models  to  represent  highly  diver¬ 
sified  actual  communities.  In  this  way  the  epidemiological  status  of  these  com¬ 
munities  could  be  tested  in  advance  of  any  real  outbreak  by  performing  appro¬ 
priate  simulations.  Alternatively,  the  models  might  be  used  to  assess  actual 
outbreaks  of  serious  disease  by  facilitating  the  estimation  of  biologically  ini- 


SIMULATION  OF  STOCHASTIC  EPIDEMICS 


257 


portant  parameters  such  as  contact  rates,  removal  rates,  and  so  forth.  From 
here  one  might  hope  to  be  able  to  construct  the  future  development  of  the 
outbreak  in  probability  terms,  indicating  a  range  of  possible  consequences  if 
various  alternative  steps  were  taken  leading  to  a  reduction  of  contact  rates  or  an 
increase  of  removal  rates,  for  example.  The  results  of  such  calculations  might 
well  indicate  what  kind  of  public  health  intervention  would  be  most  effective. 

These  ideas  are  of  course  highly  speculative,  but  indicate  some  of  the  direc¬ 
tions  in  which  further  research  might  proceed. 
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